CoCalc Logo Icon
DocsShareSupport News Sign UpSign In
Views: 100
Image: ubuntu2004
Embed | Download | Raw |
Kernel: Python 3 (ipykernel)
import random
import subprocess
import shlex
from IPython.core.display import SVG
import pyomo.environ as pyo
from pysat.solvers import Solver
from pysat.formula import CNF 
import py_svg_combinatorics as psc
from ipywidgets import widgets, HBox
from collections import Counter
from pprint import pprint
from random import randint
import numpy as np
from IPython.display import IFrame
import IPython
from copy import copy
import os
from pathlib import Path
import time

nbname = ''
try:
    nbname = __vsc_ipynb_file__
except:
    if 'COCALC_JUPYTER_FILENAME' in os.environ:
        nbname = os.environ['COCALC_JUPYTER_FILENAME']    
title_ = Path(nbname).stem.replace('-', '_').title()    
IFrame(f'https://discopal.ispras.ru/index.php?title=Hardprob/{title_}&useskin=cleanmonobook', width=1280, height=700)

Генератор случайных данных

def get_random_knapsack(n):
    B = np.random.randint(n*3, n*10)
    # предметы не влезающие в рюкзак сразу не берем
    s = np.random.randint(1, B, size=(n))
    v = np.random.randint(1, 10*B, size=(n))
    return s, v, B

s, v, B = get_random_knapsack(10)
s, v, B
(array([31, 19, 29, 22, 33, 28, 28, 43, 7, 27]), array([ 75, 120, 306, 466, 92, 268, 153, 104, 362, 197]), 47)
def get_model(s, v, B):
    m = pyo.ConcreteModel()

    m.n = s.shape[0]

    # на всякий случай, возьмем с собой
    m.s = s
    m.v = v
    m.B = B

    m.J = range(m.n)

    m.x = pyo.Var(m.J, domain=pyo.Binary)
    m.obj = pyo.Objective(expr = sum( m.v[j] * m.x[j] for j in m.J), sense=pyo.maximize)

    @m.Constraint()
    def влезает_в_рюкзак(m):
        return sum(s[j] * m.x[j] for j in m.J) <= B

    return m

m = get_model(s, v, B)
m.pprint()
1 Set Declarations x_index : Size=1, Index=None, Ordered=Insertion Key : Dimen : Domain : Size : Members None : 1 : Any : 10 : {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} 1 Var Declarations x : Size=10, Index=x_index Key : Lower : Value : Upper : Fixed : Stale : Domain 0 : 0 : None : 1 : False : True : Binary 1 : 0 : None : 1 : False : True : Binary 2 : 0 : None : 1 : False : True : Binary 3 : 0 : None : 1 : False : True : Binary 4 : 0 : None : 1 : False : True : Binary 5 : 0 : None : 1 : False : True : Binary 6 : 0 : None : 1 : False : True : Binary 7 : 0 : None : 1 : False : True : Binary 8 : 0 : None : 1 : False : True : Binary 9 : 0 : None : 1 : False : True : Binary 1 Objective Declarations obj : Size=1, Index=None, Active=True Key : Active : Sense : Expression None : True : maximize : 75*x[0] + 120*x[1] + 306*x[2] + 466*x[3] + 92*x[4] + 268*x[5] + 153*x[6] + 104*x[7] + 362*x[8] + 197*x[9] 1 Constraint Declarations влезает_в_рюкзак : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : -Inf : 31*x[0] + 19*x[1] + 29*x[2] + 22*x[3] + 33*x[4] + 28*x[5] + 28*x[6] + 43*x[7] + 7*x[8] + 27*x[9] : 47.0 : True 4 Declarations: x_index x obj влезает_в_рюкзак 
solver = pyo.SolverFactory('cbc')
solver.solve(m).write()
res = m.obj()
m.x.extract_values()
# ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: 828.0 Upper bound: 828.0 Number of objectives: 1 Number of constraints: 1 Number of variables: 10 Number of binary variables: 10 Number of integer variables: 10 Number of nonzeros: 10 Sense: maximize # ---------------------------------------------------------- # Solver Information # ---------------------------------------------------------- Solver: - Status: ok User time: -1.0 System time: 0.07 Wallclock time: 0.06 Termination condition: optimal Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available. Statistics: Branch and bound: Number of bounded subproblems: 0 Number of created subproblems: 0 Black box: Number of iterations: 0 Error rc: 0 Time: 0.17979764938354492 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0
{0: 0.0, 1: 0.0, 2: 0.0, 3: 1.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 0.0, 8: 1.0, 9: 0.0}

Приводим к задаче о разрешимости

Задача разрешима для заданного VV, если ∑u∈U′v(u)⩾V\displaystyle\sum_{u\in U'}v(u) \geqslant V.

def is_solvable(s, v, B, V, solver='cbc'):
    m = get_model(s,v,B)

    try:
        solution = pyo.SolverFactory(solver).solve(m)

        if str(solution['Solver'][0]['Termination condition']) == 'optimal' and m.obj() >= V:
            return True, m
        else:
            return False, m
    except:
        return False, m

is_solvable(s,v,B, res)[0], is_solvable(s,v,B, res + 1) [0]
(True, False)

Доказательство, что NP−comleteNP-comlete

Сведем 3SAT к нашей задаче

Сведение будем проводить в 2 этапа:

  1. Сведение 3SAT к Subset Sum Problem(SSP).

  2. Сведение SSP к Knapsack.

Формулировка SSP:

Дан массив целых чисел ai,i∈A=1..n‾a_i, i\in A = \overline{1..n}, данно целое число SS. Возможно ли выбрать подмножество A′⊆A:∑i∈A′ai=SA' \subseteq A: \displaystyle\sum_{i \in A'} a_i = S?

Сначала простая часть (пункт (2)): Построим по SSP([a],S)SSP([a], S) -> Knapsack([s],[v],B,V)Knapsack([s],[v],B,V):

Пусть B=S=VB = S = V, ∀i∈A→si=vi=ai\forall i\in A \rightarrow s_i = v_i = a_i. Тогда в задаче о рюкзаке нам нужно выбрать A′⊆A:∑i∈A′si⩽BA' \subseteq A: \displaystyle\sum_{i \in A'} s_i \leqslant B и ∑i∈A′vi⩾V\displaystyle\sum_{i \in A'} v_i \geqslant V, то есть V⩽∑i∈A′ai⩽B⟺∑i∈A′ai=S V \leqslant \displaystyle\sum_{i \in A'} a_i \leqslant B \Longleftrightarrow \displaystyle\sum_{i \in A'} a_i = S

Эквивалентность данной SSP([a],S)SSP([a], S) и полученной Knapsack([s],[v],B,V)Knapsack([s],[v],B,V) из равносильности выше.

Сведение 3SAT к SSP: Пункт (1) смотрите здесь

def cnf3_to_ssp(cnf3):
    clauses = cnf3.clauses

    var_names = set()
    for clause in clauses:
        for var in clause:
            var_names.add(abs(var))

    # Число переменных n
    n = len(var_names)
    # Число конъюкций m
    m = len(clauses)

    table = [[0 for _ in range(n + m)] for _ in range(2*n + 2*m)]
    for line in range(0, 2*n + 2*m, 2):
        table[line][line // 2], table[line + 1][line // 2] = 1, 1

    for (column, clause) in enumerate(clauses):
        for var in clause:
            line = 2 * (abs(var) - 1) if var > 0 else 2 * (abs(var) - 1) + 1
            table[line][n + column] = 1

    str_numbers = [''.join(map(str, line_of_digits)) for line_of_digits in table]
    sum_str = '1' * n + '3' * m
    return str_numbers, sum_str

def gen_random_3cnf(conjunction_count, variables_count):
    return CNF(from_clauses=psc.rand3cnf(conjunction_count, variables_count))

random.seed(1)
cnf = gen_random_3cnf(4, 3)
cnf.clauses
[(-2, -3, -1), (-3, 1, -2), (-3, 1, 2), (-3, 1, -2)]
cnf3_to_SSP(cnf)
(['1000111', '1001000', '0100010', '0101101', '0010000', '0011111', '0001000', '0001000', '0000100', '0000100', '0000010', '0000010', '0000001', '0000001'], '1113333')

В сведение мы считали, что числа выше записаны в десятиричной системе. Но так получаются слишком большие числа и pyomo умирает (pyomo работает с 32-битными интами, а у нас длина чисел n+mn+m).

Как можно немного облегчить себе жизнь: несложно доказать, что, если бы эта таблица была записана не в 10-ричной, а в 4-ричной системе, то сведение все еще будет справедливо(*). Тогда, переходя к knapsack, можно перевести эти числа в 10-ричную систему и они станут гораздо меньше (длина уменьшится в log⁡104\log_{10}4 раз). Как будет видно далее, с большими формулами pyomo все еще не справляется, но так оно начинает справляться хотя бы с небольшими.

StasFomin: Ну, если экономить, то нам надо

  • «ровно бит» на «выбор переменная или отрицание» — и переполнение там невозможно, ждем единицу,

  • и по «два бита на скобку» — где мы ждем «тройку», и для тройки тоже переполнение невозможно. — если младшие два бита переполнятся максимум мы в них получим «1», и по индукции, никакие биты для скобок не переполнятся.

  • и вот это 2*m+n — вроде в 64 битный UINT можно попробовать упихать. 24 скобки от 16 переменных, уже на что-то похоже например, вполне можно поиграть с рандомом.

Доказательство (*): Если таблица (и финальная сумма) написаны в 4-ричной системы счисления, то сведение все еще работает. От чего оно вообще может не работать? От того, что у нас в правой половине таблицы (где в каждом столбце по 5 единиц) может произойти переполнение и в другой столбец при сложении может прийти единичка. На самом деле, это ничего нам не ломает: рассмотрим сколько единиц должно быть "выбрано" для самого правого столбца: их должно быть ровно 3 (мы не сможем набрать 7 единиц, чтобы произошло переполнение и в нашем столбце была тройка: 134=71013_4 = 7_{10}, а у нас максимум 5 единиц). Значит на следующий столбец не пойдет единица. Далее доказываем по индукции, что для каждого столбца все еще нужно взять ровно 3 единицы, то есть основа сведения не поменялась.



def transform_4ns_to_10ns(initial):
    # Идем справа налево
    init = reversed(initial)

    res = 0
    current_pow_of_4 = 1

    for digit in init:
        res += int(digit) * current_pow_of_4
        current_pow_of_4 *= 4

    return res

transform_4ns_to_10ns('01001'), transform_4ns_to_10ns('123'), transform_4ns_to_10ns('100')
(65, 27, 16)

Сводим полученую SSP задачку в Knapsack (плюс перевод в десятичную)

def cnf3_to_knapsack(cnf3):
    str_numbers, str_sum = cnf3_to_ssp(cnf3)

    a, S = np.array([transform_4ns_to_10ns(num) for num in str_numbers]), transform_4ns_to_10ns(str_sum)
    return a,a, S,S

cnf = CNF(from_clauses=[(1, 2, 3), (2, 3, 4), (1, -2, -4)])
cnf3_to_ssp(cnf), cnf3_to_knapsack(cnf)
((['1000101', '1000000', '0100110', '0100001', '0010110', '0010000', '0001010', '0001001', '0000100', '0000100', '0000010', '0000010', '0000001', '0000001'], '1111333'), (array([4113, 4096, 1044, 1025, 276, 256, 68, 65, 16, 16, 4, 4, 1, 1]), array([4113, 4096, 1044, 1025, 276, 256, 68, 65, 16, 16, 4, 4, 1, 1]), 5503, 5503))

Проверяем, что наше сведение работает

Для этого генерируем случайную 3-КНФ формулу и решаем ее двумя способами: (1) используя SAT-solver и (2) превращая в экземляр нашей задачи, после чего решая Knapsack при помощи pyomo.

Проверка разрешимости 3-КНФ через преобразование к Knapsack и решение через pyomo-ЦЛП

def is_solvable_by_transforming(cnf3):
    (s,v,B,V) = cnf3_to_knapsack(cnf3)
    return is_solvable(s,v,B,V)[0]

Проверка разрешимости 3-КНФ через pysat Solver

def is_solvable_by_sat_solver(cnf3):
    from pysat.solvers import Solver

    solver = Solver(bootstrap_with=cnf3)
    return solver.solve()

Проверим что это вообще работает на тривиальных 3-КНФ.

true_cnf = CNF(from_clauses=[(1,2,3), (1,-2,-3), (1,2,-3)])
false_cnf = CNF(from_clauses=[(1,2,3), (1,2,-3), (1,-2, 3), (1,-2,-3), (-1,2,3), (-1,2,-3), (-1,-2, 3), (-1,-2,-3),])

assert(is_solvable_by_transforming(true_cnf))
assert(not is_solvable_by_transforming(false_cnf))
assert(is_solvable_by_sat_solver(true_cnf))
assert(not is_solvable_by_sat_solver(false_cnf))

Тест, который генерирует test_lengthtest\_length случайных 3-КНФ формул длины до max_conj_countmax\_conj\_count и проверяет, что наш ILP и pysat солвер выдают один и тот же ответ на задачу разрешимости


def test_many(test_length, max_conj_count, var_count=-1):
    start = time.time()
    wrong_answers = 0
    for _ in range(test_length):
        conj_count = randint(1, max_conj_count)
        if var_count == -1:
            var_count = randint(3, max(3, conj_count // 3 + 2))
        formula = gen_random_3cnf(conjunction_count=conj_count, variables_count=var_count)

        if is_solvable_by_transforming(formula) != is_solvable_by_sat_solver(formula):
            wrong_answers += 1

    return wrong_answers, time.time() - start

Теперь просто гоняем тесты

Много тестов на маленьких формулах:

test_many(500, 10)
(1, 128.98304438591003)
test_many(75, 15)
(0, 20.39103889465332)

Наше сведение хорошо решает тесты, можно предположить, что сведение легальное.

К сожалению, на больших формулах начинают происходить ошибки.

Предположение: судя по документации, Pyomo не работает с long. А числа у нас генерятся размера (p+n)∗log⁡104(p+n) * \log_{10}4, то есть при больших формулах в вычислениях pyomo может происходить переполнение

StasFomin: 'numpy.uint64' — недостаточно? Если при сведении начинаются какие-то экспоненциальные длины чисел, значит не очень это сведение-то полиномиальное. Надо посмотреть внимательно.



test_many(10, 100)
WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined>
(8, 2.047656297683716)

О многом говорит "termination condition: infeasible". Это значит, что полученная задача Knapsack неразрешима, то есть min⁡i(si)>B\displaystyle\min_{i}(s_i) > B, что явно не так (B=B = '1' * k + 3 * 'p', что больше любого sis_i). То есть произошло переполнение и BB оказалось меньше, чем на самом деле.

С этим ничего не поделать: нам не научить pyomo работать с long.



def cnf32knapsack_stas_V01(cnf3):
    n = cnf3.nv
    m = len(cnf3.clauses)

    if not( 2 * m + n <= 64):
        print(m, n, 2 * m + n)
        assert(False)

    V = int( '1' * n + '11' * m, 2)

    all_items_count = 2 * m + 2 * n
    v = np.zeros(all_items_count, dtype=np.uint64)
    C = []
    for x in range(1, n+1):
        variable_part = '0' * (x-1) + '1' + '0'*(n-x)
        for i, lit in enumerate([x, -x]):
            clause_part = ''
            for clause in cnf3.clauses:
                if lit in clause:
                    clause_part += '01'
                else:    
                    clause_part += '00'    
            row_num = 2 * (x-1) + i
            v[row_num] = int( variable_part + clause_part, 2)            

    # добиваем «мусором», чтобы выполнялось.        
    for i, clause in enumerate(cnf3.clauses):
        clause_part = '00' * i + '01' + '00' * (m - i - 1)
        int_garbage = int(clause_part, 2)
        for j in range(2):
            row_num = 2 * n + 2 * i + j
            v[row_num] = int_garbage

    s = v
    B = V

    return v, s, B, V

def print_bin(v, V):
    for i in range(v.shape[0]):
        print(f'{v[i]:#064b}')
    print(f'{V:#064b}')


cnf3 = CNF(from_clauses=[(1, 2, 3), (2, 3, 4), (1, -2, -4)])
v, s, B, V = cnf32knapsack_stas_V01(cnf3)
print_bin(v, V)
0b00000000000000000000000000000000000000000000000000001000010001 0b00000000000000000000000000000000000000000000000000001000000000 0b00000000000000000000000000000000000000000000000000000100010100 0b00000000000000000000000000000000000000000000000000000100000001 0b00000000000000000000000000000000000000000000000000000010010100 0b00000000000000000000000000000000000000000000000000000010000000 0b00000000000000000000000000000000000000000000000000000001000100 0b00000000000000000000000000000000000000000000000000000001000001 0b00000000000000000000000000000000000000000000000000000000010000 0b00000000000000000000000000000000000000000000000000000000010000 0b00000000000000000000000000000000000000000000000000000000000100 0b00000000000000000000000000000000000000000000000000000000000100 0b00000000000000000000000000000000000000000000000000000000000001 0b00000000000000000000000000000000000000000000000000000000000001 0b00000000000000000000000000000000000000000000000000001111111111 
B, V
(1023, 1023)
ok, m = is_solvable(s,v,B,V)
ok, m
(True, <pyomo.core.base.PyomoModel.ConcreteModel at 0x7f92e68798b0>)
m.x.extract_values()
{0: 0.0, 1: 1.0, 2: 0.0, 3: 1.0, 4: 1.0, 5: 0.0, 6: 1.0, 7: 0.0, 8: 1.0, 9: 1.0, 10: 1.0, 11: 0.0, 12: 1.0, 13: 1.0}
def is_solvable_anyhow(m, V):
    try:
        OK = False
        m.max_value = pyo.Constraint(rule = m.obj >= V)
        for solver in ['cbc', 'scip']:
            solution = pyo.SolverFactory('cbc').solve(m)
            if str(solution['Solver'][0]['Termination condition']) == 'optimal':   #and m.obj() >= V
                return True, m
            print(f'Pyomo called {solver} cannot solve it :(')    

        if not OK:
            file_mps = 'sat-knapsack.mps'
            m.write(file_mps)
            scmd = f'scip -c "read {file_mps} optimize quit"'
            output = subprocess.check_output(shlex.split(scmd)).decode('utf-8')
            if 'optimal solution found' in output:
                print(f'Directly called SCIP solve it :)')    
                return True, m
            return False, m    
    except Exception as ex_:
        print(ex_)
        return False, m

def is_feasible_SAT(cnf3):
    solver = Solver(bootstrap_with=cnf3)
    return solver.solve()

def generate_cnf3():
    m = randint(1, 29)
    n = randint(3, 62 - 2 * m)
    if not( 2 * m + n <= 64):
        print(m, n, 2 * m + n)
        assert(False)

    cnf3 = CNF(from_clauses=psc.rand3cnf(m, n))
    assert(len(cnf3.clauses) == m)
    return cnf3

def run_tests(t):
    tests_count_positive = tests_count_negative = 0
    while tests_count_positive < t or tests_count_negative < t:
        cnf3 = generate_cnf3()
        is_feasible_SAT_ = is_feasible_SAT(cnf3)
        if is_feasible_SAT_:
            if tests_count_positive >= t:
                continue
            tests_count_positive += 1
        else:    
            continue
            if tests_count_negative >= t:
                continue
            tests_count_negative += 1

        print(f'Testing on CNF3(n={cnf3.nv}, m={len(cnf3.clauses)}), feasible={is_feasible_SAT_}')
        v, s, B, V = cnf32knapsack_stas_V01(cnf3)
        m = get_model(s, v, B)
        solvable, m = is_solvable_anyhow(m, V)
        if is_feasible_SAT_ != solvable:
            print(cnf3)
            assert(False)
    
    print("OK")

run_tests(10)
Testing on CNF3(n=7, m=9), feasible=True Testing on CNF3(n=13, m=18), feasible=True Testing on CNF3(n=8, m=27), feasible=True WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> Pyomo called cbc cannot solve it :( WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> Pyomo called scip cannot solve it :( 
def wtf_cbc_01():
    v = np.array([2251799814012928, 2251804108652816, 1125905275551744,
                  1125917371924480,  562949953421312,  562949953422400,
                    281474976714756,  281474976710657,  140737488355328,
                    140737572241408,   70368744177664,   70385924050965,
                    35184372154368,   35185450025984,   17592454479888,
                    17592202821953,    8796093022208,    8796093088768,
                      4399388688384,    4398113619972,    2199023517696,
                      2199023255552,    1099511644160,    1103877898240,
                      549757140992,     549755813888,     292058828864,
                      274877906944,     137438953472,     137438953728,
                        68724736000,      68719476736,      17179869184,
                        17179869184,       4294967296,       4294967296,
                        1073741824,       1073741824,        268435456,
                          268435456,         67108864,         67108864,
                          16777216,         16777216,          4194304,
                            4194304,          1048576,          1048576,
                            262144,           262144,            65536,
                              65536,            16384,            16384,
                              4096,             4096,             1024,
                              1024,              256,              256,
                                64,               64,               16,
                                16,                4,                4,
                                  1,                1], dtype=np.uint64)
    V = 4503599627370495 
    s = v
    B = V
    cnf3 = CNF(from_string="p cnf 16 18\n-2 -6 14 0\n2 -12 -1 0\n-7 10 2 0\n10 -2 8 0\n-10 -5 -12 0\n-8 -5 -2 0\n16 -12 -7 0\n16 14 13 0\n13 1 11 0\n-9 1 7 0\n16 13 12 0\n4 -6 14 0\n-3 -9 -7 0\n-1 -15 -8 0\n-8 -3 14 0\n-6 -1 8 0\n-6 -10 4 0\n-4 -8 -6 0\n")                     
    solver = Solver(bootstrap_with=cnf3)
    is_feasible_SAT = solver.solve()
    #sat_solution = solver.get_model()
    sat_solution = [-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, 13, -14, -15, -16]
    solvable, m = is_solvable(s, v, B, V)

    m = get_model(s, v, B)
    solution1 = pyo.SolverFactory('cbc').solve(m)
    print(solution1.write())
    solution2 = pyo.SolverFactory('scip').solve(m)
    print(solution2.write())
    solution2
    ...

wtf01()
WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> WARNING: Loading a SolverResults object with a warning status into model.name="unknown"; - termination condition: infeasible - message from solver: <undefined> # ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: inf Upper bound: None Number of objectives: 1 Number of constraints: 0 Number of variables: 0 Number of binary variables: 68 Number of integer variables: 68 Number of nonzeros: 0 Sense: maximize # ---------------------------------------------------------- # Solver Information # ---------------------------------------------------------- Solver: - Status: warning User time: -1.0 System time: 0.0 Wallclock time: 0.0 Termination condition: infeasible Termination message: Model was proven to be infeasible. Statistics: Branch and bound: Number of bounded subproblems: None Number of created subproblems: None Error rc: 0 Time: 0.11304640769958496 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0 None # ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Lower bound: -inf Upper bound: 4503599627370490.0 Number of objectives: 1 Number of constraints: 0 Number of variables: 68 Sense: unknown # ---------------------------------------------------------- # Solver Information # ---------------------------------------------------------- Solver: - Status: ok Message: optimal solution found Termination condition: optimal Id: 0 Error rc: 0 Time: 1.86 Gap: 0.0 Primal bound: 4503599627370490.0 Dual bound: 4503599627370490.0 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 1 number of solutions displayed: 1 - Status: unknown Objective: No values Variable: No values Constraint: No values None 
def wtf_scip_01():
    cnf3 = CNF(from_string="p cnf 21 11\n7 -1 6 0\n16 -13 -9 0\n8 -19 -7 0\n-12 -21 9 0\n-9 10 -14 0\n5 -11 -13 0\n3 -8 20 0\n-15 19 4 0\n10 -21 20 0\n19 17 2 0\n-3 -18 17 0\n")
    solver = Solver(bootstrap_with=cnf3)
    is_feasible_SAT = solver.solve()
    sat_solution = solver.get_model()
    #sat_solution = [-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, 13, -14, -15, -16]
    v, s, B, V = cnf32knapsack_stas_V01(cnf3)
    solvable, m = is_solvable(s, v, B, V, solver='cbc')

    m = get_model(s, v, B)
    solution2 = pyo.SolverFactory('scip').solve(m)
    # m.hint_constraints = pyo.ConstraintList()
    # for lit in sat_solution:
    #     xn = abs(lit)
    #     if lit > 0:
    #         m.hint_constraints.add(expr =  m.x[ 2*(lit-1) ] == 1)
    #         m.hint_constraints.add(expr =  m.x[ 2*(lit-1)+1 ] == 0)
    #     else:
    #         m.hint_constraints.add(expr =  m.x[ 2*(lit-1) ] == 0)
    #         m.hint_constraints.add(expr =  m.x[ 2*(lit-1)+1 ] == 1)
    solution3 = pyo.SolverFactory('scip').solve(m)
    ...

wtf_scip_01()
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) /var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb Cell 59 line 1 ----> <a href='vscode-notebook-cell://xn--17-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb#Y151sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a> wtf_scip_01() /var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb Cell 59 line 1 <a href='vscode-notebook-cell://xn--17-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb#Y151sdnNjb2RlLXJlbW90ZQ%3D%3D?line=16'>17</a> m.hint_constraints.add(expr = m.x[ 2*(lit-1)+1 ] == 0) <a href='vscode-notebook-cell://xn--17-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb#Y151sdnNjb2RlLXJlbW90ZQ%3D%3D?line=17'>18</a> else: ---> <a href='vscode-notebook-cell://xn--17-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb#Y151sdnNjb2RlLXJlbW90ZQ%3D%3D?line=18'>19</a> m.hint_constraints.add(expr = m.x[ 2*(lit-1) ] == 0) <a href='vscode-notebook-cell://xn--17-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb#Y151sdnNjb2RlLXJlbW90ZQ%3D%3D?line=19'>20</a> m.hint_constraints.add(expr = m.x[ 2*(lit-1)+1 ] == 1) <a href='vscode-notebook-cell://xn--17-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/maximum-knapsack.ipynb#Y151sdnNjb2RlLXJlbW90ZQ%3D%3D?line=20'>21</a> solution3 = pyo.SolverFactory('scip').solve(m) File ~/hard-problems-formulations/.venv/lib64/python3.11/site-packages/pyomo/core/base/var.py:1033, in IndexedVar.__getitem__(self, args) 1031 def __getitem__(self, args): 1032 try: -> 1033 return super().__getitem__(args) 1034 except RuntimeError: 1035 tmp = args if args.__class__ is tuple else (args,) File ~/hard-problems-formulations/.venv/lib64/python3.11/site-packages/pyomo/core/base/indexed_component.py:646, in IndexedComponent.__getitem__(self, index) 644 if isinstance(index, EXPR.GetItemExpression): 645 return index --> 646 validated_index = self._validate_index(index) 647 if validated_index is not index: 648 index = validated_index File ~/hard-problems-formulations/.venv/lib64/python3.11/site-packages/pyomo/core/base/indexed_component.py:874, in IndexedComponent._validate_index(self, idx) 867 raise KeyError( 868 "Cannot treat the scalar component '%s' " 869 "as an indexed component" % (self.name,) 870 ) 871 # 872 # Raise an exception 873 # --> 874 raise KeyError( 875 "Index '%s' is not valid for indexed component '%s'" % (idx, self.name) 876 ) KeyError: "Index '-4' is not valid for indexed component 'x'"